#################################################################
## Tyler Pratt ##################################################
## "Deference and Hierarchy in International Regime Complexes" ##
## Code for bootstrapped substantive effects ####################
#################################################################

## Bootstrap Table 3, model 1 to gauge substantive effects (shown in Fig. 4)

boot.results <- as.data.frame(matrix(nrow = 1000, ncol = 5))
colnames(boot.results) <- c("majors_diff", "IO2_weight", "IO2_tech", 
                             "Binding_Tech", "Decision")
eqn <- "glm(deference_prop ~ majors_diff + IO2_weight + 
                IO2_Technical + Binding_Technical + DecisionMaking_Ease_Dir +
                regime + t + I(t^2) + I(t^3),
            data = new.data, family = binomial(link = 'logit'))"

Def <- read.csv("Deference_Replication.csv")
Dir.Dyads <- as.character(unique(Def$Dir.DyadID))

set.seed(9999)
for(i in 1:1000){
  new.dyads <- sample(Dir.Dyads, length(Dir.Dyads), replace = TRUE) # clustered sampling by directed dyad
  new.data <- Def[0,]
  for(j in 1:length(new.dyads)){
    new.data <- rbind(new.data, Def[which(Def$Dir.DyadID==new.dyads[j]),])
  }
  result <- eval(parse(text = eqn))
  
  dat1 <- dat2 <- new.data
  dat1$majors_diff <- 0
  dat2$majors_diff <- 3
  boot.results$majors_diff[i] <- mean(predict(result, newdata=dat2, type="response") - 
                                         predict(result, newdata=dat1, type="response"), na.rm=T)
  dat1 <- dat2 <- new.data
  dat1$IO2_weight <- 0
  dat2$IO2_weight <- 1
  boot.results$IO2_weight[i] <- mean(predict(result, newdata=dat2, type="response") - 
                                        predict(result, newdata=dat1, type="response"), na.rm=T)
  dat1 <- dat2 <- new.data
  dat1$IO2_Technical <- 0
  dat2$IO2_Technical <- 1
  boot.results$IO2_tech[i] <- mean(predict(result, newdata=dat2, type="response") - 
                                      predict(result, newdata=dat1, type="response"), na.rm=T)
  dat1 <- dat2 <- new.data
  dat1$Binding_Technical <- 0
  dat2$Binding_Technical <- 1
  boot.results$Binding_Tech[i] <- mean(predict(result, newdata=dat2, type="response") - 
                                          predict(result, newdata=dat1, type="response"), na.rm=T)
  dat1 <- dat2 <- new.data
  dat1$DecisionMaking_Ease_Dir <- 0
  dat2$DecisionMaking_Ease_Dir <- 1
  boot.results$Decision[i] <- mean(predict(result, newdata=dat2, type="response") - 
                                      predict(result, newdata=dat1, type="response"), na.rm=T)
  print(i)
}
write.csv(boot.results, "bootBL_Def2017.csv")

plot(colMeans(boot.results)[c(5,4,3,2,1)], 1:5, bty="n", pty="sq",
     yaxt = "n", ylim = c(0.4, 5.6), ylab="",
     xlab = "Change in Deference Intensity", pch = 16, cex = 1,
     xlim = c(-0.045, 0.11),
     col = c(rep("blue", 3), rep("red", 2)))
abline(v=0, lty=2)
lines(x=quantile(boot.results$Decision, probs = c(0.025, 0.975)), y=c(1,1), col = "blue", lwd=2)
lines(x=quantile(boot.results$IO2_weight, probs = c(0.025, 0.975)), y=c(2,2), col="blue", lwd=2)
lines(x=quantile(boot.results$IO2_tech, probs = c(0.025, 0.975)), y=c(3,3), col = "blue", lwd=2)
lines(x=quantile(boot.results$IO2_weight, probs = c(0.025, 0.975)), y=c(4,4), col="red", lwd=2)
lines(x=quantile(boot.results$majors_diff, probs = c(0.025, 0.975)), y=c(5,5), col = "red", lwd=2)
text(x=colMeans(boot.results)[5], y=1.15, "Decision-Making Difference", col = "blue", cex=0.9)
text(x=colMeans(boot.results)[4], y=2.15, "Binding-Technical Pair", col = "blue", cex=0.9)
text(x=colMeans(boot.results)[3], y=3.15, "Technical IO", col = "blue", cex=0.9)
text(x=colMeans(boot.results)[2], y=4.15, "Weighted Voting", col = "red", cex=0.9)
text(x=colMeans(boot.results)[1], y=5.15, "Major Power Difference", col = "red", cex=0.9)

